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ABSTRACT 

There is ample observational evidence that the star formation rate (SFR) surface density, Esfr, is 
closely correlated with the surface density of molecular hydrogen, Eh 2 • This empirical relation holds 
both for galaxy- wide averages and for individual > kpc sized patches of the interstellar medium (ISM) , 
but appears to degrade substantially at a sub-kpc scale. Identifying the physical mechanisms that 
determine the scale-dependent properties of the observed Eh 2 — Esfr, relation remains a challenge from 
a theoretical perspective. To address this question, we analyze the slope and scatter of the Eh 2 — ^SFR 
relation using a set of cosmological, galaxy formation simulations with a peak resolution of ~ 100 pc. 
These simulations include a chemical network for molecular hydrogen, a model for the CO emission, 
and a simple, stochastic prescription for star formation that operates on ~ 100 pc scales. Specifically, 
star formation is modeled as a Poisson process in which the average SFR is directly proportional to the 
present mass of H2 . The predictions of our numerical model are in good agreement with the observed 
Kennicutt-Schmidt and Eh 2 — Esfr relations. We show that observations based on CO emission 
are ill suited to reliably measure the slope of the latter relation at low (< 20 M Q pc -2 ) H2 surface 
densities on sub-kpc scales. Our models also predict that the inferred Eh 2 — Esfr relation steepens 
at high H2 surface densities as a result of the surface density dependence of the CO/H2 conversion 
factor. Finally, we show that on sub-kpc scales most of the scatter in the relation is a consequence 
of discreteness effects in the star formation process. In contrast, variations of the CO/H2 conversion 
factor are responsible for most of the scatter measured on super-kpc scales. 

Subject headings: galaxies: evolution - galaxies: formation - starsdormation - methods: numerical 



1. INTRODUCTION 

The formation of stars is far from being a well 
understood and solved problem. The complex interplay 
between star formation, the formation of molecular 
clouds, and the physical processes operating within 
the interstellar medium (ISM), such as turbulence, 
magnetic fields, self-gravity, or feedback from the 
stellar population, poses many formidable challenges 
to the proper interpretatio n and modeling of th e 
star formation process (e.g., IMcKee fc OstTiken 12001 . 
Fortunately, observations have revealed a number of 
empirical relations which provide some guidance for the 
development of a theoretical model. This includes the 
well-studied Kennicutt-Schmidt relation, a correlation 
betwee n the su r face d ensity of the ne u tral I SM and 
E SFR (|Schmidtl [i95l IKennicuttl [l989l floM) . The 
realization that star formation is more clo sely correlated 
with molecular gas than ne utral gas (iTalbotl 119801 : 
IWong fc Blitl l200l but cf. IKennicuttl 11989) led to 
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Bigiel et al.l 120111: iFeldmann et all 1201 H iSchruba et all 
20111: iRahman et al.ll201lL 120121 ). 

The properties of the Eh 2 — Esfr. relation are tied 
to the small scale (scales of molecular clouds and be- 
low) connection between star formation and the sup- 
ply in form of molecular gas. Therefore, measuring the 
slope and scattei0 of the Eh 2 — Esfr relation allows to 
probe potentially relevant physical mechanisms that are 
involved in the star formation process. This requires, of 
course, that biases in the observational tracers used to 
derive Eh 2 and Esfr are properly accounted for. The 
main goal of this paper is therefore to study and iso- 
late the importance of both various physical mechanisms 
and observational biases for the slope and scatter of the 
Eh 2 — Esfr relation. 

Our theoretical predictions are based on high 
resolution, cosmological galaxy formation simula- 
tions equipped with subgrid models to estimate 
H2 abundances and SFRs. We study observational 
biases by post-processing our simulations with a 
model for the J = 1 — > line emission of CO 
(|Feldmann et al.l l2012t paper I from now on). CO 
line emission is the most comm only used tra cer to 
measure H 2 abundanc e s (e.g., iWilson et al.l 
Scoville fc Sanders! 119871: iBrown fc Vanden Boutl 1 1991 
Young et al.lll995t iReean et al.ll2001l: iHelfer et aLl l2003 
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5 This is a somewhat imprecise, but common, terminology. Here 
(and in the rest of the paper) "slope" and "scatter" of the Ejj 2 — 
Esfr relation refer, more precisely, to the slope and to the standard 
deviation of the residuals (estimated from a linear regression) of the 
log 10 Sh 2 - log 10 Esfr relation. 
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llvison et al.1 120111) . Detecting H 2 m emission is 
difBcult due to the high energy gap between the ground 
state and the lowest excited levels of H2 and the lack of 
a permanent electric dipole moment of the H2 molecule. 
Our CO emission model allows us to determine the 
CO/H2 conversion factor (also called the X- factor or 
Xco) which is formally defined as 



X 



CO 



T^co 



(1) 



i.e., as the ratio between H2 column density, Nn 2 , and the 
CO velocity integrated intensity, Wco, of the J = 1 — > 
rotational transition. As shown in paper I, the X-factor 
depends on the conditions of the ISM, in particular on the 
dust-to-gas ratio and the gas surface density. Hence, a 
crucial question that we address in this paper is whether 
the use of a constant value for Xco, as done in many 
observational studies, introduces significant biases in the 
estimates of the slope and scatter of the £h 2 — ^sfr rela- 
tion. Similarly, we use stellar population synthesis mod- 
eling to compute the FUV and Ha fluxes of our model 
galaxies in order to assess whether the use of these star 
formation tracers leads to observational biases. 

From a theoretical perspective, the slope of the Eh 2 — 
SgFR relation provides us with a means to distinguish 
between (and potentially rule out) different models of 
star formation. In general, the slopes of the Eh 2 — ^sfr 
relation and that of the underlying relation between the 
spatial densi ties of SFRs and H 2 (jGuibert et al.l [19781 : 
lTaiboti[i980[) . the pn 2 — p* relation, coul d be different 
due to variations in the gas scale height ( |Talbotl Il971f) 
or ch anges in the density distribution (jFeldmann et al] 
1201 If) . However, if the P h-> _ P* relation is linear , so is 
the Sho — ^sfr relat ion (jSchave fc Dalla Vecchial [20081 : 
IFeldmann etal|[20rl . 

We note that a linear relation has been interpreted 
as evidence that star formation occurs in dense clumps 
within molecular clouds and that the number of such 
clumps scales with the total amount of molecular gas 
(jBigiel et al.l [20 08). An alternative possibility that has 
been suggested is that the timescale for star formation 
is controlled by small scale processes, e.g., stellar feed- 
back, ambipolar diffusion, etc. which do not vary much 
from one location to the next, and hence ensure a steady 
(when averaged over sufficiently large sp atial and tem- 
poral scales) conversion of H2 into stars (Wong & Blitz 
2002). In contrast, a slope of ~ 1.3 can be explained 
by models that are based on a constant star formation 
efficiency per free-fal l time ()Krumholz fc McKcej 120051 : 
iKrumholz et al.ll2009f ). Star formation driven by cloud- 
cloud collisions could lead to even steeper slopes (|Taskerl 

limb. 

While the slope of the £h 2 — Ssfr relation has been 
the subject of a number of observational studies, there 
is, unfortunately, no clear consensus yet on its ex- 
act value. While many studies find slopes that are 
close to unity (~ 0.8 - 1.3, e.g.. IWong fc Blitzl 1 20021: 
Komugi et al.|[20ql i lBTgiel et al.ll2008HBlanc et aljbooa 
Bigiel et all 120111: ISchruba et all 120111: IRahman et al l 
20111 I2012I) . a number of work s prefer a steeper slope , 



e.g., ~ 1.4 dHever et all [200l IKennicutt et all 
1.6 (|Thilker et al.ll2007l) . ~ 1.2-1.9 (ILiu et al.l 
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2011) 



eluding choices in the fitting methodology, how star for- 
mation maps are corrected for dust extinction and diffuse 
emission, the resolution of the observations, or the range 
in surface densities over which the fit is done. In partic- 
ular, assumptions about the amount of diffuse emission 
present in the star formation maps and its subsequent 
treatment in the data analysis haye a s t rong impact on 
the derived slope (IKennicutt et aLll2007l; ILiu et al.ll2011t 



o pe 

IRahman et al.ll2011f) . However, if the fitting is restricted 
to regions of sufficiently high Eh 2 , then the impact of any 
diffuse emission is minimized and the £h-> — Ssfr re lation 
is found to be close to linear (Rahman et al. 2012). 

Besides the slope, the normalization of the Eh 2 — 
Ssfr relation is an important parameter that needs 
to be matched by the theoretical models. Fortu- 
nately, studies that find an approximatively linear 
slope of the £h 2 — Ssfr relation also find little vari- 
ation of the normalization (£h 2 /£sfr ~ few Gyr) 
with galaxy mass or other global galaxy properties 



s or other global galaxy propertiet 
2001; Wong & Blitz 2002'; Bigiel et al 



2010; lBigiel et al.l l2011: Bo latto et al 



Various factors may contribute to these differences, in- 



(e.g.-lBraine et aL 
| 2008t iGenzel et al 
2011). A notable exception is the claim that the nor- 
malization depends on the specific star formation rate 
(jSaintonge et al.l 120111 ). although there is a potential 
worry that this result is, to some extent at least, an arti- 
fact of the large scatter in the Eh 7 ~ £-sfr relation an d 
the use of correlated quantities ( Rahman et al.l [2012). 
Driven by these observational findings we adopt in this 
paper a (stochastic) star formation model that is based 
on a linear pu 2 — P* relation with a constant gas depletion 
time of a few Gyr. 

The scatter in the £h 2 — ^sfr relation has received sig- 
nificantly less attention compared to the slope or the nor- 
malization. The reason is probably that the scatter de- 
pends even more on the specific details of the respective 
observational survey and the data analysis procedures. 
The finding that the scatter increases with increasing 
spatial resolution of a survey is sometimes used to argue 
that star form ation scaling relations "break down" on 
small scales (IMomose et all 120101 : IQnodera et al.l 120101 : 
ISchruba et atfe oiO). However, we will argue in this pa- 
per that the proper interpretation is that on small scales 
the intrinsically stochastic nature of the scaling relations 
becomes simply more evident. 

A first step to understand and quanti fy the different 
contri butors to the scatter was made by IFeldmann et all 
(poTlfh Here, we extend our previous analysis in several 
ways. We consider the scatter that results from spatial 
X-factor fluctuations. We also analyze the scatter re- 
lated to star formation tracers (Ha and FUV). Finally, 
we highlight the importance of stochasticity in the star 
formation process and propose a simple model of star for- 
mation that allows a classification of the various mecha- 
nisms that contribute to the scatter in the Sh 2 — £<sfr 
relation. 

The outline of the paper is as follows. In we present 
the detail of our numerical approach, including the set-up 
of the simulations ( §2.1|) . the star formation model fi j2.2|) . 
details of sub-grid physics ( §2.3[) . and the modeling of 
the H2 and star formation tracers ( §2.4j) . We then show 
in ij3.ll that our simulations are able to reproduce the 
observed Kennicutt- Schmidt and Sh 2 — Ssfr relations. 
Next, in Q3.2[ we discuss observational claims of a non- 
linear pn 2 — P* relation based on slope measurements of 
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TABLE 1 

Details of the numerical simulations 



label 


box size 


Ax 


m DM 




[Sim, S2a, fib, h, os\ 


fixed ISM 


comments 


MW-fid 


6 com. Mpc h _1 


65 pc 


9 x 10 5 M Q hr 


i 


[0.3, 0.7, 0.043, 0.7, 0.9] 


yes 


fiducial SF model 


MW-dt6 


6 com. Mpc h 


65 pc 


9 x 10 5 Mq hr 


l 


[0.3, 0.7, 0.043, 0.7, 0.9] 


yes 


Ai SF = 10 6 yr 


MW-dt5 


6 com. Mpc h 


65 pc 


9 x 10 5 Mq hr 


l 


[0.3, 0.7, 0.043, 0.7, 0.9] 


yes 


At SF = 10 5 yr 


MW-sl2 


6 com. Mpc h" 1 


65 pc 


9 x 10 5 Mq hr 


l 


[0.3, 0.7, 0.043, 0.7, 0.9] 


yes 


slope-2 SF model 


HZ-csm 


25.6 com. Mpc h — 1 


97 pc 


1 X 10 6 Mq h" 


l 


[0.28, 0.72, 0.046, 0.7, 0.82] 


no 


down to z = 1.8 


HZ-fid 


25.6 com. Mpc h" 1 


97 pc 


1 x 10 6 Mq hr 


l 


[0.28, 0.72, 0.046, 0.7, 0.82] 


yes 


fiducial SF model 



Note. — The first five columns provide the name, box size, numerical resolution, and the adopted cosmological parameters of each 
zoom-in simulation. The penultimate column states whether the simulation is run with gas-to-dust ratios and interstellar radiation fields 
fixed to the corresponding values in the s olar neighborhood. The parameter Aigp in the last column refers to the average time between 
individual star formation events, see i|2,2l The fiducial SF model uses Atgp = 10 yr. 



the Eh 2 — Ssfr relation at low £h 2 ■ Then, in §373] we 
analyze the importance of the X-factor for high redshift 
galaxies with high gas surface densities. We address the 
origin of the scatter in the £h 2 — Ssfr relation in §3.41 
In 5?] we discuss our star formation model in the context 
of various observational constraints. Finally, in iJSJ we 
summarize our results and conclusions. 

2. METHODOLOGY 
2.1. Simulations 

All simulations have been run with the Eulerian hy- 
drody namics -I- N-body code ART (jKravtsov et al.lll997l 
2002) that uses an adaptive mesh refinement (AMR) 
technique to increase the resolution selectively in speci- 
fied regions of interest. We also use the standard method 
of embedding these regions in layers of subsequently 
lower dark matter resolution to further reduce the com- 
putational cost, but s till capture the impact of large scale 
tidal fields correctly (|Katzl)l99lHBertschingerll200lD . 

Simulation MW-fid focusses its computational re- 
sources on a Lagrangian region that encloses five virial 
radii of a MW-sized halo at z — (total mass ~ 10 12 
Mq) in a 6 Mpc h -1 box. The simulation is started from 
cosmological initial conditions with the parameters given 
in Table [T] It is run fully self-consistently down to red- 
shift z = 4. At this point "fixed ISM conditions" are 
imposed on the simulation as described in §2.31 and it is 
continued for additional 600 Myr before it is analyzed. 
By then the high-resolution Lagrangian region harbors a 
large disk galaxy with a virial mass of^4.2xl0 11 Mq 
and several less massive galaxies. More details on the 
setup of the MW-fid simu l ation can be found els ewhere 
(|Gnedin fc Kravtsovl [201lt iFeldmann et alJl20TH paper 
I). 

The simulations MW-dt5, MW-dt6, and MW-sl2 are 
spawned from the z = 4 snapshot of MW-fid and con- 
tinued for additional 600 Myr with fixed ISM conditions. 
The first two simulations are run with reduces values 
of A^sf, the average time scale between individual star 
formation events, see ^2.21 but are otherwise identical to 
MW-fid. The run MW-sl2 uses a modified star formation 
model, see £13.21 for details. 

Simulation HZ-csm refines on seven regions within a 25 
Mpc h^ 1 box. By z = 0, these regions have collapsed to 
halos in the 10 11 — 10 13 Mq mass range. The simulation is 
started at z = 100 and is continued down to z = 1.8 fully 
self-consistently. The s etup of the HZ-cs m simulation 
is discussed in detail in iZemp et al.l (|2012[ ). Simulation 
HZ-fid is spawned from the z = 2 snapshot of HZ-csm. 



It is continued with fixed ISM conditions for additional 
200 Myr (down to z = 1.8) to allow the gas to react to 
the changes in the ISM properties and to reach its new 
equilibrium state. 

The main properties of our simulations are summarized 
in Tabled] 

2.2. Star formation model 

Since individual resolution elements in our simulations 
correspond to (at best) GMC scales, it is clear that we 
are not able to follow in any realistic detail the formation 
and evolution of individual bound star clusters, let along 
individual stars. Hence, our approach is to marginalize 
over the complexities involved in star formation by de- 
scribing star formation on a statistical level. Such an ap- 
proach will have its limitations, of course, and we are not 
aiming at reproducing observations perfectly, but rather 
try to isolate and explain some of the trends in the ob- 
servational data. 

There are several lines of evidence suggesting that a 
statistical description of star formation is a plausible 
ansatz, at least on the level of star clusters and in rela- 
tively quiescently star-forming galaxies. One is the obser- 
vation that the scaling between Ssfr and £h 2 lS roughly 
linear and holds over a wide range of surface densities 
and g alaxy metallicities (|Bigiel et al.ll2008t iGenzel et al.l 
2010), although it may break down in strongly out- 
of-equilibrium situat i ons such a major gala xy mergers 
(jTevssier et al.1 120101 iBournaud eFall 1201 lh . Another 
one is provided by the similarity of the shape of the mass 
function of young star clusters among non-starbursting 
galaxies ()Porte gics Zw art et al.l I2010I) . which indicates 
that the conversion of interstellar gas into star clusters 
proceeds on average in the same way in different galaxies 
with different global properties. Also, the existence of a 
correlation between the total number of clusters and the 
luminosity of the most luminous cluster in a galaxy can 
be understood i n a purely st atistical formation scenario 
of star clusters (|Larse nl l2002Th 

Our mo del assumes that H 2 mass is a good tracer of 
the SFR dPelupessy et al.l 120061: iRobertson fc Kravtsov 



2008 1 : IGnedin et al.1 l2009TlPapadopoulos fc Pelupessyl 
201(1 IGnedin fc Kravtsovl^Ollt IFeldmann et al.l 1201 ll: 



Kuhlen et al.l|201 2l) and can be summarized by the fol- 
lowing three statements: (1) star formation is a stochas- 
tic process and the number of individual star formation 
eventf@ in a given time interval is described by a ho- 
mogeneous Poisson process, (2) the average stellar mass 

6 The proper interpretation of an individual star formation event 
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formed per unit time is proportional to the present H2 
mass, and (3) the factor of proportionality, 7"^ , is a con- 
stant. Here, average refers to the ensemble average over 
independent patches of the ISM of the same size and with 
the same H2 masses, which may, however, differ in other 
properties. In other words, we marginalize over most of 
the GMC parameters which may be potentially relevant 
for star formation and keep only the explicit dependence 
of the SFR on H2 mass. 

We give a detailed description of the model in the ap- 
pendix. Here, we only summarize the properties that we 
will need later in the paper. 

The discreteness of star formation leads to scatter in 
the Sh 2 — Ssfr relation because the SFR that is real- 
ized in a given region and over a given period differs 
from the ensemble average SFR tied to the present H2 
mass. The scatter increases (up to a point) with in- 
creasing average time between individual star formation 
events, AisF, but decreases with increasing life-time of a 
given star formation tracer At*. In fact, the scatter de- 
pends on these quantities only in form of the dimension- 
less ratio At»/A£gF, the average number of individual 
star formation events during At* . 

Our model contains three parameters: the average H2 
depletion time scale, Td ep , the spatial scale on which our 
model operates, Isf, and AtsF- We fix Td cp to 2.9 Gyr 
in order to fit the observed normalization of the Esfr - 
Eh 2 relation, see £13.11 We stress that Td op is not the H2 
gas depletion time of individual star forming molecular 
clouds, but an average depletion time that includes the 
non star forming molecular gas. Our model is assumed 
to operate on scales of Isf = 60 — 100 pc scales, which 
is the peak resolution of our numerical simulations, see 
Table III Our default value for At SF is 10 Myr. This 
choice is based on the following considerations. 

First, the crossing time (or the free-fall time) of a 
molecular region should be a lower limit on the time pe- 
riod between two bursts of star formation in that region. 
This time i C ross ~ L/a can be estimated from the Larson 
relation between cloud size L and cloud v elocity disper- 
sion a (Larson 198lt [Solomon et all 119871) . It scales as 
~ 1 Myr (L/pc) - 5 and, hence, t cross ~ 8 - 10 Myr for re- 
gions of size Isf = 60 — 100 pc. Since the largest molecu- 
lar cloud complexes in the Milky Way have radii of ~ 100 
pc ([Dame et al.lll986|) . this line of reasoning cannot be 
used to constrain AisF on larger spatial scales. Despite 
having a similar numerical value on ~ 100 pc scales, Atgp 
and Across are quite different physical quantities. In fact, 
Across, a measure of the duration of an individual star for- 
mation event, increases with scale, while A^sf, the pe- 
riod between successive star formation events, decreases 
with scale due to the smoothing effect of spatial averag- 
ing on Poisson noise. 

Secondly, we have checked that with the choice A^sf = 
10 Myr individual star formation events in the simula- 
tions cover a broad range of masses up to 10 6 M©. The 
embedded star clusters that form during a star formation 
event of total mass 10 6 M Q will have masses up to, but 
not exceeding, ~ 10 5 M Q for a reasonable choice of the 
slope of the embedded cluster mass function (~ 2.3—2.4). 
This compares well with the fact that, at least in the 

is as the simultaneous formation of a number of embedded star 
clusters in a given region in space (~ 60—100 pc in our simulations). 



Milky Way, young massive clusters exceeding 10 5 M & 
are rare ([Portegies Zwart et al.|[2010h . 

Finally, AtsF ~ 10 Myr is consistent with the observed 
relation b etween SFR and max imum embedded star clus- 
ter mass ([Wcidncr et al. 2004]). We note that these con- 
siderations do not rule out a moderately larger value, e.g, 
Ai S F ~ 20 Myr. 

2.3. Further sub-grid modeling 

All simulations include a photo-chemical network, 
metal enrichment from supernova (type la and type II), 
but no thermal energy injection, optically thin radiative 
cooling by hydrogen (including H2), helium, and metal 
lines, and 3D radiative transfer of UV radi ation. The de- 
tails o f the implementation can be foun d in lGnedin et al.l 
(|2009h and lGnedin fc Kravtsov! (pOll . Here, we give a 
brief recount. 

The photo-chemical network in ART follows the for- 
mation and destruction of the five major atomic and 
ionic species of hydrogen and helium. The formation of 
H2 on dust grains and the destruction of H2 via photo- 
dissociation in the Lyman- Werner bands are taken into 
account. The transfer of ionizing and non-ionizing UV 
radiation from s tellar sources is compu ted in the OTVET 
approximation ([Gnedin fc Abell feOOL). Radiative stellar 
feedback is important for the heating and cooling bal- 
ance of the gas and for the abundance of H2 and CO. 
Unlike the re-ionized intergalactic medium, the dense in- 
terstellar medium within galaxies may well be opaque to 
ionizing photons of all but the nearest stars and radiative 
transfer effects cannot be neglected. 

Most of ou r simulations are run with "fixed ISM con- 
ditions" (see lGnedin fc Kravtsov! [20TlT ) . First of all this 
means that the dust-to-gas ratio Duw , which is normally 
assumed to scale linearly with the local gas metallicity, 
is kept fixed at a value that corresponds to3 Z@, and 
we will denote this as Z?mw = 1- The gas-to-dust ra- 
tio is a crucial parameter that enters the formation rates 
and the dust shielding of molecular hydrogen and also 
our CO emission model. The metallicities of the self- 
consistently enriched gas are still used to compute the 
cooling rates that enter the hydrodynamical solver. This 
is done in order to avoid numerical artifacts such as sud- 
den increases in the gas accretion rates resulting from 
excess, non-equilibrium cooling. 

Furthermore, in the fixed ISM runs, the normaliza- 
tion of the radiation field at 1000 A is fixed to Jmw 
= 10 6 photons cm -2 s _1 sr _1 eV _1 , a val ue typical for 
the solar neighborh ood in the Milky Way (|Drainel[l978l ; 
IMathis et all 11983). Wc use the notation Umw = L 
where Umw is the intensity of the radiation field at 1000 
A in units of Jmw- We stress that only the normaliza- 
tion of the radiation field computed with the OTVET 
solver is fixed. The shape of the radiation spectrum is 
not modified. 

2.4. Postprocessing: CO, FUV and Ha emission 

7 In this paper Zq refers to the metallicity of the solar neighbor- 
hood. Specifically, Z Q = 0.02 or 12 + log 10 (O/H) = 8.92, which 
is somewhat larg er than the metallicity of t he Sun according to 
recen t estimates HAllende Prieto et all 120011 ; lAsplund et ail 120031 , 
[20091) . 
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Fig. 1. — Relation between the molecular gas (neutral gas) surface density Sjj 2 (£hi+h 2 ) an d the surface density of the star formation 
rate EgpR as predicted by our fiducial SF and CO emission model for Milky- Way like ISM properties (Dmw = 1, ^mw = !)■ The left 
panel shows the Eu 2 vs Esfr relation, the right panel the classical Kennicutt-Schmidt relation (£hi+h 2 vs ^SFr)- The spatial resolution 
is chosen to match as closely as possible the resolution of the observational studies (Bigicl ct al. 2008, 2011). The gas surface densities 
include a factor 1.36 that accounts for the presence of Helium. The red dot-dashed line shows the median EgFR (time averaged over the 
past 20 Myr) for a given "true" H2 or total gas surface density (as computed in the simulation). The blue solid line shows instead the 
median of Ssfr as a function of the "inferred" H2 or total gas surface density, i.e., a gas density in which Su 2 is inferred from the CO 
intensity, as predicted by our model, using the galactic X-factor. The magenta hashed region indicates the typical scatter in Ssfr (25 and 
75-th percentiles in the left panel, 16 and 84-th percentiles in the right panel, ). The black dashed line (best fit) an d the contour lines 
(containing 90% and 50% of the data) in the left panel are observational measurement of the Eh 2 - Ssfr relation by [Bigicl ct all l|201U ) 
on kpc scales. Th e contour line in the right panel shows the distribution of sub-kpc sized patches in the sample of nearby galaxies by 
IBigiel et af] ||2008T) , see their Fig. 7. The figure demonstrates that our modeling of the star formation and the CO emission is consistent 
with the observed relations between Ssfr and the molecular and neutral gas surface density, respectively. 



The J = 1 — > 12 C 16 emission is computed as de- 
scribed in paper I. In short, a sub-grid model for the CO 
emission is constructed based on the results of a suite 
of small scale magneto-hy drodynamical ISM simulations 
(Gl over fc Mac Lowll2011l ). This model contains two free 
parameters. 

One is the CO brightness temperature that is related 
to both the temperature of the CO emitting gas and 
the temperature of the cosmic microwave background 
(CMB). In i|3J] $LI and flwe study the £h 2 - Ssfr 
in the local Universe and we adopt a gas temperature 
of 10 K (a typical temperature of molecular clouds in 
the Milky Way). The corresponding brightness temper- 
ature is 6.65 K. In ^3.31 we predict the CO emission 
for galaxies at z ~ 2. We assume that the increase 
in the CMB temperature at those redshifts is compen- 
sated by an increase in the gas temperature from 10 K 
to 14.5 K, such that the brightness temperature remains 
approximatively constant. Such a moderate increase in 
gas temperature is consistent with the detailed modeling 
of the gas temperature in non-starbursting high redshift 
galaxies using photon-dissociation regions c odes coupled 
with a semi-analytic galaxy evolution model (jLagos et al.1 
[20T1 . 

The other free parameter is the scaling of the CO line 
width (either a constant line width or a virial scaling). 
By default we show results for the virial line width scal- 
ing, which is probably the more realistic of the two (see 
Fig. 3 of paper I), but if relevant we will point out the 



changes that result from the assumption of a constant 
line width. The sub-grid model is applied to the high- 
est refined resolution elements in the simulation. The 
contributions from these individual, ~ 60 — 100 pc sized 
resolution elements are then combined in the optically 
thin limit to derive the CO emission from larger regions. 

The FUV and Ha emission of each stellar partic le is 
computed with Starburst-99 ([Leitherer et al.lT2010[ ) as- 
suming solar metallicity, high mass loss Geneva tracks. 
We checked that switching to, e.g., Padua tracks does 
not affect any of our results in a significant way. The 
Ha luminosities are directly taken from the Starburst- 
99 output, while the broad-band FUV luminosities are 
computed using the Galex FUV transmission curve 
(|Morrissev et al.l [20051 ) and the UV spectra provided by 
Starburst-99. 

3. RESULTS 

3.1. Actual and inferred star formation relations 

It has been shown dGnedin et al . 2009; Krumhol z et al.l 
2009; IGnedin fc Kravtsovi 12010) that an H 2 -based star 
formation prescription is able to reproduce the relation 
between the surface densities of neutral gas and SFRs 
(the Kennicutt-Schmidt relation) both for galaxies in 
the local Universe and for galaxies at high redshift. In 
fact, the connection between star formation and molecu- 
lar hydrogen provides a simple physical interpretation 
for the drop in the SFR surface densities at low gas 
surface densities (Robertso n fc Kravtsovi |20"08() . In this 
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Fig. 2. — Same as the left panel in Fig[T] but for a SF model in which p t a p| . The left panel uses surface densities measured on kpc 
scales, while the right panel shows the corresponding results for spatial averaging scales of 250 pc. This SF model leads to a steepening 
of the inferred £jj 2 — Ssfr relation at sufficiently high H2 surface densities (> 20 Mq pc~ 2 ). However, this steepening is suppressed for 
Sjj 2 < 20 Mq pc - 2 , especially when measured at sub-kpc spatial resolution, as a result of the increase of Xqo with decreasing Eh 2 at 
low surface densities (paper I). It will therefore be challenging to extract an unbiased estimate for the slope of the relation between SFR 
and gas density from CO observations of MW-like galaxies at low surface densities. 



picture the drop is a manifestation of a transition be- 
tween a neutral and a molecular hydrogen phase. Specif- 
ically, above a characteristic gas surface density the gas 
is shielded from the interstellar radiation field by a suffi- 
ciently large dust optical depth and the molecular phase 
prevails. The characteristic surface density depends on 
the dust-to-gas ratio, the strength of the interstellar ra- 
diat ion field, and the density structure of the ISM, see, 
e.g. lGnedin fc Kravtsov! (pOifl ). 

However, this interpretation of the Kennicutt-Schmidt 
relation neglects a potentially relevant detail. While the- 
oretical models can predict H2 surface densities directly, 
observations often rely on CO observations to infer H2 
surface densities. A crucial test of our understanding of 
the Kennicutt-Schmidt relation is therefore whether this 
agreement still holds even if we take the effects of the 
CO/H2 conversion factor into account. In other words, 
we need to compare theory and observations on an equal 
footing, i.e., using H2 surface densities that are derived 
from CO emission in both cases. 

Our star formation model contains the gas depletion 
time Tdep as a free parameter. This parameter is the con- 
version factor between H2 mass and the ensemble average 
SFR, see i )2.2[ but, since we assume Td ep =const, it also 
corresponds to the normalization of the £h 2 — ^sfr rela- 
tion. We find that a depletion time Td cp of 2.9 Gyr leads 
to an inferred (CO based) depletion time of ~ 2.35 Gyr 
(including Helium, Bigi el et aLll2011[ ) and, thus, to excel- 
lent agreement between predictions and observations, see 
Fig. QJ,. We note that the proper choice of Td ep is largely 
degenerate with the value of the X- factor. Specifically, 
our CO emission model prefers a median Xco value for 
a Milky Way like ISM that that is ~ 65% larger than 
the value X q o.mw = 2 x 10 20 cm' 2 K _1 km -1 s used by 
iBigiel et al.l (|2011[ ). explaining the need for a somewhat 



larger H 2 depletion time. 

Fig. QJ, shows the inferred (CO based) £h 2 — ^sfr re- 
lation based on our simulation MW-fid. It is an approx- 
imatively linear relation over two orders of magnitude in 
H2 surface density. This is not an entirely obvious result 
since, at least on sufficiently small scales, the CO/H2 
conversion factor is strongly dependent on surface den- 
sity (paper I). However, it turns out that the spatial av- 
eraging over a large set of regions with different Xqq 
values erases most of this surface density dependence on 
kpc and larger scales. Hence, the actual (linea r) slo pe of 
the Eh 2 — Ssfr relation is recovered (but see ^3.31) . 

We show our predictions for the Kennicutt-Schmidt 
relation in Fig. [T]d, finding good agreement between our 
theoretical predictions and observations of the inferred 
(CO based) Kennicutt-Schmidt relation. We thus con- 
clude that a star formation model that depends linearly 
on the abundance of H2 is consistent with both the 
observed Kennicutt-Schmidt relation and the observed 
Eh 2 — ^sfr relation of normal star-forming galaxies in 
the local Universe. 



3.2. 



J SFR 



relation 



The slope of the Sh 2 

The slope of the Sh 2 — Ssfr relation has been a subject 
of significant debate over the last years. While a number 
of studies find a slope of about unity, other observational 
works favor steeper slopes. It has been argued that these 
discrepancies arise from observational obstacles such as 
diffuse emission in the infrared, or uncertainties in the 
dust absorption that make it challenging to obtain reli- 
able SFR estimates, especially in regions of low star for- 
mation activity and gas surface d ensity (|Liu et alj|2011t 
iRahman et~aTTl2011fc iLerov et al.ll2012l) . 

A varying CO/H2 conversion factor is another of those 
potential obstacles, but has been largely neglected since 
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there is no general observational handle on this quantity 
in different galaxies. The aim of this section is to test the 
implications for the observed Eh 2 — Esfr relation if the 
p* — Pn 2 relation is highly non-linear and the X-factor 
is taken into account. In the following we assume that 
any observational uncertainties related to the estimates 
of SFR (e.g., diffuse infrared emission, dust absorption, 
etc.) are under control and corrected for. We focus in- 
stead on the role of the CO/H 2 conversion factor. 

For this test we use simulation MW-sl2, which uses the 
same general approach to star formation as presented in 
§2.2[ but with a density dependent gas depletion time. 
Specifically, we assume that Td cp oc 1//9h 2 which implies 
that (on average) p* oc . The factor of proportionality 
is chosen such that the gas depletion time, as estimated 
from the inferred (CO based) Eh 2 — Esfr relation mea- 
sured on kpc scales, is ~ 2.3 Gyr at Eh 2 = 10 Mq pc~ 2 . 

The results of this test are shown in Fig. [2] While 
the actual Eh 2 — Esfr relation is very steep with a 
slope n ~ 2, the inferred (based on CO observations) 
Eh 2 — Esfr relation flattens significantly at Eh 2 % 20 
Mq pc~ 2 . This effect is already pronounced on kpc 
scales, but becomes very strong on scales of 250 pc. In 
other words, we predict that the H2 gas depletion time 
derived from CO emission would appear to decrease, not 
increase, with decreasing surface density at Eh 2 < 20 
Mq pc -2 on sufficiently small scales. This effect is a con- 
sequence of the anti-correlation between Xco and Eh 2 at 
low gas surface densities, caused by t he presence of large 
amounts of CO-dark molecular gas (iWolfire et all 120101 : 
Krumholz et al.l 12011b INaravanan et al.l l2011bt paper I; 
Shettv et al.ll2011fc iNaravanan et al.ll2012t) . 

For 20M Q pc" 2 < E H2 < lOOM pc" 2 , on the other 
hand, the X-factor plays only a small role and the steep 
slope (n ~ 2) is recovered. At such gas surface densities 
the diffuse emission in the star fo rmation maps is als o 
expected to be less of an obstacle (|Rahman et al.| [2012). 
Hence, the good news is that a robust measurement of 
the slope of the Eh 2 — Esfr relation based on CO data 
will be possible in areas of moderately high gas and SFR 
surface density. At low Ejj 2 , however, and especially on 
sub-kpc scales, an accurate determination of the slope on 
the basis of CO observations will be difficult. 

3.3. Comparison with z ~ 2 galaxies 

Studies of the E H2 — Esfr relation that are based on 
samples of galaxies from both the local Universe and from 
high redshift indicate that the slope is close to (but not 
quite) linear over a many orders of magnitude of Eh 2 
(|Genzel et al.l 12010). In particular, if one selects only 
non-interacting galaxies, the slope of the Eh 2 — Esfr re- 
lation is approximatively ~ 1.1 — 1.2. This, however, 
appears to be driven by a steepening at high surface 
densities, because the slope measured on i ntermediate 
H2 su rface densities is much closer to unity (jBigiel et al.l 
[20081) . 

In our analysis we focus on non-interacting galaxies, 
that presumably form stars in an approximatively steady 
state. We do not include interacting or merging galaxies 
in our discussion, as su ch objects can be scattered off 
the Eh ^ — Esfr relation (|Genzel et al.ll2010llDaddi et al.l 
I2010bfl . 

While there are several ways to explain a non-linear 



slope of the Eh 2 — Esfr relation, we demonstrate below 
that our star formation and CO emission model coupled 
with an actually linear Eh 2 — Esfr relation does predict a 
slightly non-linear slope of the inferred (CO-based) Eh 2 — 
Esfr relation at high Eh 2 , consistent with observations. 
In essence, we argue that the observed super-linear slope 
can be understood as an X-factor effect. 

We note that the following predictions rest on the as- 
sumption that there is no significant trend of the CO 
line ratios or the CO brightness temperature with gas 
surface density. Our predictions are specifically for the 
,7 = 1—5-0 rotational transition line of CO, while galaxies 
at z ~ 1 — 2 are typically observed in CO line emission 
resulting from J = 2->lorJ=3->2 transitions. 
Systematic trends in the emission line ratios could thus 
modify the slope of the Eh 2 — Esfr relation. A system- 
atic variation of the brightness temperature with surface 
density would have a similar effect. So far there is no 
clear observational evidence that either of these assump- 
tions is violated in steady-state, normally star forming 
galaxies at z ~ 1 — 2. 

Our predictions for the inferred (CO based) Eh 2 —Esfr 
relation at z ~ 2 are shown in Fig. [3^. Galaxies with 
inferred molecular gas densities Eh 2 > 100 M Q pc~ 2 are 
shifted off the actually linear Eh 2 — Esfr relation. The 
origin of this offset is the small, but systematic, variation 
of the X-factor with surface density. 

As discussed in paper I, at high Ejj 2 the CO emission 
from a small ISM patch (~ 20 — 100 pc) ceases to scale 
linearly with the gas column density (and thus the H 2 
surface density) due to the increased optical thickness 
at the line center of the CO emission line. This effect 
is somewhat, but not fully, compensated by an increase 
in the width of the emission line, so that overall there 
is a remaining increase of the X-factor with increasing 
H2 surface density. For this reason, the use of a constant 
CO/H2 conversion factor leads to systematic shifts in the 
inferred Eh 2 — Esfr relation. At intermediate H2 surface 
densities (the precise range depends on the spatial scale) 
the conversion factor is essentially constant and, hence, 
in this case the measured s lope is predicted to be very 
close to linear, as observed (|Bigiel et al.H2008l) . 

Star forming galaxies at higher redshifts typically have 
lower Z than galaxies of a s i milar mass in the local 
Universe, e.g.. iMaiolino et al.l (|2008D . Our fully self- 
consistent cosmological simulation HZ-csm predicts that 
metallicities are only ~ 0.5^0 at z ~ 2. Hence, Dmw ~ 
0.5 and one may wonder whether the combination of a 
high redshift, low metallicity sample and a low redshift, 
high metallicity sample could explain the observed super- 
linear relation. In order to test the importance of the 
metallicity dependence of CO/H2 conversion factor we 
rerun the HZ-csm simulation with fixed Milky- Way like 
ISM conditions (simulation HZ-fid, see N2.ip . 

Fig. [5b shows that galaxies lie along the observed, 
slightly super-linear Eh 2 — Esfr relation, even if they 
had Dmw — 1- We therefore conclude that the super- 
linear slope is caused primarily by the scaling of Xco 
with H2 surface density, and that metallicity and dust- 
to-gas ratio variations have only a small effect. 

We stress again that our CO emission model rests on 
a number of assumptions as pointed out above. If either 
of these assumptions were broken, an alternative expla- 
nation for the super-linearity of the Eh 2 — Esfr relation 
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Fig. 3. — The Sjj 2 — Ssfr relation as predicted by a fully cosmological, hydrodynamical simulation at z ~ 2 (simulations HZ-csm, 
HZ-fid). Each panel shows the predictions of the "inferred" £r 2 — ^SFR relation, i.e., the one where H2 surface densities are derived from 
12 CO emission maps using a galactic conversion factor -Xqo.mw = 2 x 10 20 K — 1 cm -2 km -1 s. (Left) the interstellar radiation fields 
and metallicities as computed self-consistently within the simulation at 2 = 1.8. The H2 weighted mean metallicity at this epoch is only 
~ 0.5 Zq but spans a broad range from (1-u) 0.25 to 0.75 Zq . This results in a value of Xco that is above the galactic conversion factor and 
hence shifts the median of the inferred Sjj 2 — Ssfr relation slightly toward lower surface densities when compared with the observations 
by Bigicl ct al. (2011) for galaxies in the local Universe. The more interesting result is, however, that the slope of the Epj 2 — Ssfr relation 
appears to steepen, in particular at high EgpR. (Right) The simulation is restarted at a slightly earlier epoch and continued for ~ 200 
Myr down to z = 1.8, but this time with dust-to-gas ratios and UV radiation fields fixed to D^w = 1 and Umw = b respectively. The 
steepening of the Ejj 2 — Ssfr relation remains visible and is therefore not a result of changes in the dust-to-gas ratios or interstellar 
radiation fields. Symbols and lines are as in Fig. [Tp. In addition, the gray dashed line shows the fit to the observed Sjj 2 — SgpR relation 
based on a large sample of low and high-z galaxies by Genzcl ct al. (2010). The red triangles mark the individual positions of simulated 
galaxies with stellar masses exceeding 10 10 Mq. All simulation predictions are based on our fiducial SF and CO model. The latter assumes 
a virial scaling of the CO line width. The red arrow at the top indicates the median shift in the inferred H2 column density of the simulated 
galaxies if the CO line width would be fixed to a constant value of 3 km s —1 . The figure shows that ga laxies with h i gh gas or SFR 
surface densities appear to deviate from a linear Ejj 2 — SgpR relation, consistent with the observations of lGenzel et al.l 1 120101). despite 
the fact that the underlying relation between SFR and H2 mass is perfectly linear. The super-linear slope (~ 1.1 — 1.2) of the inferred 
Sjj 2 — SgpR relation is caused primarily by the increase of Xqq with increasing Sjj 2 at high gas column densities. Our results apply 
only to galaxies that are in an equilibrium mode of star formation, not to starbursting galaxies. In the latter environments our CO model 
becomes unr eliable and the galaxy - wide ratio between total molecular gas and the for SF relevant dense molecular gas (n > 10 4 cm -3 ) 
may change IIGao fc So lomon 200l; ILada et al.ll20H ; IPapadopoulos et al.ll2012l) . 



would be required. 

Systematic variations of the CO line ratios or the CO 
brightness temperature with gas surface density could be 
responsible for a difference between the intrinsic slope of 
the Eh 2 — Ss fr relation and the slope derive d from CO 
observations ([Narayan an et all 1201 lal |2012|) . Alterna- 
tively, a variation of the CO line ratios or brightness tem- 
perature with redshift could lead to a systematic z depen- 
dence of the normalization of the observed Sh 2 — ^sfr 
relation. This can produce an artificial trend with sur- 
face density if a galaxy sample is used in which Ejj 2 
(and Ssfr) strongly correlate with galaxy redshift. Fi- 
nally, the Er 2 — Ssfr relation may actually get steeper 
at high column densities. For instance, it has been sug- 
gested that star formation becomes more efficient at high 
column densities because external pressure on molecular 
clouds sh ifts the balance between g ravity and turbulent 
support (|Krumholz fc McKeel 120051) . Upcoming obser- 
vations with the Atacama Large Millimeter Array will 
hopefully enable us to disti nguish between our model and 
these alternatives (e.g., see IFu et alJ[20Tl . 

3.4. The scatter in the Sh 2 — SgpR relation 



There are many effects and processes that could, in 
principle, contribute to the scatter in the £h 2 — ^sfr 
relation. Clearly, scatter can arise from (1) uncertainties 
related to the method of estimating SFRs, (2) uncer- 
tainties related to the estimation of H2 masses and 
surface densities, (3) a possible non-linearity of the star 
formation process, and (4), any systematic uncertainties 
in the observables that were not accounted for. We 
will focus in this paper on the scatter sources (1) and 
(2). The poten t ial ro le of (3) is discussed in detail in 
iFeldmann et al.l ()2011l ). We do not attempt to model 
sources that fall under category (4), since they are not 
intrinsic but depend on the specifics of the observational 
survey. 

SFR estimates: The stellar mass that was formed over 
some past time interval will, in general, not coincide 
with the SFR that is expected based on the present 
H2 m ass. In other words, as shown in IFeldmann et al.l 
((20T1 . the use of a time- averaged SFR as an estimator 
of the ensemble- average SFR introduces scatter. The fol- 
lowing (not necessarily distinct) mechanisms fall under 
this category: 
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• discreteness of star formation: star formation oc- 
curs in individual star formation events, i.e., is clus- 
tered in time, 

• stochasticity of star formation: star formation re- 
lations on small scales hold only in an (ensemble) 
average sense. 

• fluctuations in the H2 abundance: H2 densities and, 
thus, the ensembl e average SFRs may fluctu ate on 
short time scales ([Glover fc Mac Lowll2007t ). while 
the SFRs derived from tracers are smoothed be- 
cause of the inherent time averaging, 

• evolutionary processes: the conversion efficiency 
from gas to stars may ch ange over the life- 
time of molecular cloud s (iMurravl 120111: cf. 
Feldmann fc Gnedin l 120111 : see also iSchruba et al.l 
20101 : lOnodera et al.ll2010f) . 

• imperfect tracers: observational tracers of SFRs 
(e.g., Ha luminosities) may provide only approxi- 
mate estimates of time-averaged SFRs because the 
proper conversion factor is not known exactly (e.g., 
it depends on the precise star formation history). 

Observations indicate that most stars form in quan- 
tized units (embedded star clusters, lLada fc Ladall2003l) 
and hence that star formation is a discretized process. 
While this inevitably makes star formation a stochas- 
tic process (in the above sense), it does not mean that 
all stochasticity in star formation arises from discrete- 
ness effects. For instance, the SFR in a molecular cloud 
depends on more than just its molecular mass. Mag- 
netic fields, cosmic ray density, or the virialization state 
of the cloud will play a role to some extent. Marginal- 
izing over these additional control parameters will give 
rise to an apparent stochasticity of star formation. In 
the formalism of H2.21 the gas depletion time Td cp be- 
comes a function of these additional parameters and its 
replacement with some appropriately averaged, constant 
depletion time leads to the appearance of stochasticity. 
For instance, a star formation efficiency that evolves over 
the lifetime of molecular clouds can be interpreted in this 
sense. Here, the age of the cloud is the additional control 
parameter that determines Td ep . 

Our numerical models account for scatter caused 
by fluctuations in the H2 abundance, the use of ob- 
servational SFR tracers, the time discreteness of star 
formation, and any incidental stochasticity, but do not 
include the potential contributions from any additional 
control parameters. 

H2 mass estimates: As mentioned above, uncertainties 
in the estimates of the H2 surface densities contribute 
to scatter in the £h 2 — Ssfr. relation. Hence, variations 
in the CO/H 2 conversion factor are a potential source of 
scatter for observational studies that are based on CO 
emission. As discussed in paper I, the X-factor can vary 
significantly, even for a fixed dust-to-gas ratio and H2 
column density, because £h 2 depends on the product of 
H2 mass fraction and total gas density, while the CO 
emission depends on the latter but not the former. 

Systematic variations of observables: We measure the 
scatter under the condition that the dust-to-gas ratio 



and the interstellar radiation field are kept fixed (in the 
sense of S}273]). This is an important point since the X- 
factor depends strongly on the former quantity (but only 
weakly on the latter, see paper I). Hence, variations in 
the dust-to-gas ratio lead to systematic modulations of 
conversion factor and, if unaccounted for, to scatter in 
the Sh 2 — Ssfr relation. 

Often mctallicity variations within a given galaxy 
are a strong function of galacto-centric radius (|Searlel 
fl97lh and can be large from one galaxy to another. 
Hence, such systematic X-factor variation will appea r 
as a galaxy-to-galaxy scatter (|Schruba et al.l 120 111 ). 
Since the amount of scatter that is created in this 
way depends on the particular selection function of 
the galaxy sample, we do not include galaxy-to-galaxy 
scatter in this analysis. Our predictions should therefore 
be compared with observations based on samples of 
galaxies with approximatively the same dust-to-gas ratio. 

In Fig. [4^, we plot the scatter of the Sh 2 — Ssfr relation 
as a function of spatial averaging scale and separate the 
contributions that result from the use of time-averaging 
tracers of star formation and uncertainties in Xcoj re- 
spectively. The scatter is computed from all regions with 
an (inferred) H2 column density between 10 and 100 
M© pc -2 , a CO velocity integrated intensity equal to or 
larger than 0.2 K km s , and a minimum SFR surface 
density of 3 x 10~ 4 Mq yr _1 kpc -2 . These limits are cho- 
sen to roughly mimic typical values encountered in obser- 
vational studies and it is clear that the exact numerical 
predictions will depend to some extent on these limits. 
In particular, the choice of the minimal Ssfr is crucial, 
since the scatter is measured in log Eh 2 — log Ssfr space 
and, if not removed, regions with very low star formation 
would contribute enormously to the scatter (even worse, 
regions with zero star formation would make the scatter 
formally infinite). 

Our numerical modeling predicts that X-factor varia- 
tions induce a scatter of the order of ~ 0.1 — 0.2 dex. 
Fig. shows that this scatter may be relevant on scales 
~ kpc and above, but on smaller scales the total scat- 
ter is primarily due to the use of time averaged SFRs, 
at least for our fiducial choice Aisp = 10 Myr. Again 
we stress that this analysis assumes that variations in 
the dust-to-gas ratio and CO brightness temperature are 
small or accounted for. Fig. 0^ also shows that the scat- 
ter that results from the use of time-averaged SFRs is a 
strong function of scale. It reaches ~ 0.5 dex at ~ 100 
pc, but is only ^0.1 dex at kpc scales. 

Furthermore, there is little difference between the use 
of FUV luminosity-based SFRs and the use of the ac- 
tual time-averaged SFR over the last 20 Myr. Hence, a 
varying FUV-to-stellar mass conversion factor (see im- 
perfect tracers above) contribu tes little to t he sc atter. 
This rules out the suggestion bv lLerov et al.l f|2012T ) that 
the variation of the FUV luminosity over the lifetime of 
a single stellar population (SSP) dominates the scatter, 
at least on scales of ~ 100 pc and above. In fact, if 
the luminosity-weighted SFRs differ little from the time- 
averaged SFRs, e.g., if SFRs are constant, then the lu- 
minosity evolution of the tracer becomes completely ir- 
relevant. This can be seen from the following simple 
analysis. 
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The total spectral luminosity from a Lagrangian vol- 
ume element at time t is given as 



L v {t) 



SFR(t')<M< - t')dt' . 



Fig. 4. — Scatter in the Sh 2 — Ssfh, relation as function of spatial averaging (resolution) scale for a galaxy with Milky-Way like ISM 
properties (Dmw = b Umw = 1). In both panels the scatter is measured over the range 10 Mq pc -2 < £jj 2 < 100 Mq pc - 2 . Regions 
with Sspr < 3 X 10~ 4 Mq yr~ 1 kpc -2 or Iqo < 0.2 K km s — 1 are excluded from the analysis. (Left) The blue dot-dashed line shows 
the scatter that arises when Eh 2 is know exactly, but SFR are inferred from FUV luminosities. The blue dotted lines marks the analogous 
result when SFR are derived from the stellar masses formed within the last 20 Myr. Scatter in Xqq at fixed CO emission leads by itself 
to a scatter in the inferred Sjj 2 — £sfr relation of the order of 0.1-0.2 dex and is shown as the magenta hashed region. The lower and 
upper boundaries of this region correspond to the cases of virial scaling of the CO line width vs constant line width, respectively (see text). 
Fluctuations in Xqo are n °t an important source of scatter on sub-kpc scales (at fixed dust-to-gas ratio and interstellar radiation field), 
but become increasingly relevant on scales of ~ kpc and above. Finally, the black horizontally hashed region shows the combined scatter 
that takes into account both the scatter in Xqq and the scatter associated with the estimations of SFRs. (Right) This panel shows how the 
scatter depends on the SF tracer (FUV, Ha, or simple time averaged SFR) and on assumptions about the stochasticity of the SF process 
(see legend). The critical parameter is the average time AtgF between SF events at a given site within the galaxy (see text). Specifically, 
the upper 4 lines show the scatter as derived from the various SF tracers for AtgF = 10 Myr (our fiducial value), while the two lines just 
below correspond to Atgp = 1 Myr. The gray hashed region shows the scatter for a run with Atgp = 0.1 Myr. This scatter results from the 
mismatch in time scales between SFRs that are averaged over the lifetime of a particular tracer (4 Myr - lower boundary; 20 Myr - upper 
boundary) and H2 masses that are observed at a given instant. The figure shows that stochastic effects play a crucial role in determining 
the overall scatter in the Ejj 2 — Ssfr relation. Furthermore, modulo Xqo effects, the scatter decreases with scale £ roughly as a power 
law oc l~ a , with a 0.5 — 0.7, consistent with the findings and interpretation given by [Feldmann ct al. (2011). 

the time discreteness parameter Atgp. We find that 
Ha-based SFR estimates lead to more scatter in the 
Eh 2 — Sspr relation than the use of FUV flux as a tracer. 
Hence, the star formation tracer with the shorter lifetime 
(Ha) leads to larger scatter. Similarly, when we estimate 
the SFR based on the actual stellar mass formed within 
the past 4 Myr and the 20 Myr, we find that the use 
of a shorter averaging time leads to more scatter in the 
Eh 2 ^ Sgppj relation. 

The averaging timescales of 4 and 20 Myr correspond 
roughly to the luminosity weighte d timescales of Ha and 
FUV emission (|Lerov et al.ll2012T) . It is therefore not en- 
tirely surprising that the scatter predictions computed 
using these time-averaged SFRs are similar to the pre- 
dictions that use Ha and FUV based tracers. It demon- 
strates that Ha and FUV based tracers can, to a good 
degree of approximation, be treated as a top-hat filter 
with a width of ~ 4 Myr and ~ 20 Myr, respectively. 
This correspondence will break, however, if the scales are 
small enough and the lifetimes of the particular tracer 
short enough such that luminosity weighted SFRs and 
time-averaged SFRs begin to differ substantially. For Ha 
based SFR estimates this appears to happen on scales of 
< 400 pc, while for FUV based tracers the effect is small 
even on spatial averaging scales of ~ 100 pc. 



Here, <fiv(t) is the spectral luminosity from a SSP of unit 
mass and age t. This can be rewritten as 

L v {t) = <SFR> „ (2) 

where 



E u = 4> v (t')dt', and 
Jo 

1 r°° 

(SFR) „ {t) = — SFR(t - t')Mt')dt' 

&v Jo 

are the total spectral energy emitted by an SSP of unit 
mass and the luminosity weighted SFR, respectively. 
Hence, if time-averaged and luminosity weighted SFRs 
trace each other closely, then the constant E v is the per- 
fect (i.e., scatter- free) conversion factor between tracer 
luminosity and time-averaged SFR. The validity of this 
statement does not depend on the form of <fi u . 

We show in Fig. |4Jd that the scatter in the £h 2 — Ssfr 
relation depends on both the star formation tracer and 
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Fig. 0J) also shows that the scatter depends on the 
average time between star formation events AisF- As 
expected a shorter Atsp means that individual star for- 
mation events involve less stellar mass but occur at a 
higher rate, which reduces the scatter. The sharp drop 
in the scatter when Aigp is reduced to 1 Myr proves that 
much of the scatter in our fiducial AigF = 10 Myr model 
arises from a single source, the discreteness of star forma- 
tion. Hence, observational estimates of the scatter can 
be used to put tight constraints on the value of A^sf- In 
addition, a systematic observational study of the scatter 
in the £h 2 — Esfr relation as function of ISM environ- 
ment would allow to determine whether (and how) Atgp 
depends on ISM properties. 

In order to assess how much of the total scatter is 
caused by time variations in the H2 abundance (and 
not related to time discreteness of star formation), we 
also show the scatter for a run with At$F = 0.1 Myr. 
This timescale is close to the smallest dynamical time 
step in our simulation and, hence, effectively eliminates 
any discreteness (beyond that dictated by the simulation 
time step) in the star formation model. We then mea- 
sure the scatter in the true Eh 2 — Esfr relation using 
time-averaged SFRs (using 4 Myr and 20 Myr as av- 
eraging times). In this way we isolate the scatter that 
is caused by the observational actuality that SFRs are 
time-averaged quantities, but H2 masses are observed at 
a particular instant. The scatter that results in this way 
is relatively small, see Fig. 2Jd. On super-kpc scales it is 
dominated by the scatter that is caused by Xqo fluctu- 
ations and on sub-kpc scales by the scatter due to the 
discreteness of star formation. However, we point out 
that time variations in the H2 abundance couple in a 
non-linear way to the Poisson noise of individual star 
formation events (see appendix). Hence, they contribute 
to the scatter caused by the discreteness of star formation 
and, hence, cannot be negl ected. 

In iFeldmann et al.l (|201lH we found that noise inserted 
by hand on small scales leads to scatter in the Ejj 2 ~£sfr 
relation that decreases roughly as power law oc l~ a , 
where a « 0.5, with increasing spatial averaging scale 
I. Fig. 3b shows that the scatter due to the stochastic 
nature of star formation follows this scaling approxima- 
tively. In particular, it is clearly less steep than a a = 1 
scaling which would be the naive expectation if the gas 
were arranged in disk o f uniform density. As discussed in 
IFeldmann et al.l (|2011[ ) the scaling deviates from a = 1 
because the density distribution of the ISM, determined 
by turbulence, is far from being uniform. 

How does the scatter that is predicted for our fiducial 
choice A^sf = 10 Myr compare with observations? Such 
a comparison is not straightforward since observational 
estimates of the scatter depend on choices in the method- 
ology. For instance, the treatment of diffuse emission 
does not only affect the inferred slope of the £h 2 — Ssfr 
relation, but also the scatter. Furthermore, the mea- 
sured scatter depends on the surface density range of the 
fit. With these caveats in mind we will now compare our 
predictions to observational studies that give quantitive 
e stimates of the s catter at a given scale. 

iRahman et al.l (|2011[) infer a scatter of about 0.3-0.4 
dex for SFRs based on Ha luminosities on ~ 0.5 kpc 
scales. The scatter is lower (~ 0.1 — 0.3 dex) when they 
use tracers with longer lifetimes (FUV + 24/xm). Our 



model predicts a scatter of ~ 0.4 dex for Ha based SFRs 
and a scatter of ~ 0.2 dex if FUV luminosities are used 
to trace star formation on ~ 0.5 kpc scales. This quan- 
titative agreement is a further justification of the choice 
AfoF = 10 My r 

iVerley et al.1 (|2010t ) investigate how the scatter in- 
creases with increasing spatial resolution. Their Fig. 4 
shows a scatter of ~ 0.4 dex on 360 pc scales, ~ 0.35 dex 
on 720 pc scales, and ~ 0.3 dex on 1.4 kpc scales. This is 
a somewhat shallower scaling than our predictions (~ 0.4 
dex on 500 pc scales, ~ 0.25 dex on kpc scales and 0.1- 
0.15 dex on super-kpc scales), see Fig. [4]b. However, their 
scatter is computed by an iterative clipping method that 
would tend to underestimate the scatter if the scatter is 
la rge, i.e., on smaller s cales . 

ISchruba et all IpOlOh compare the gas depletion times 
tco and th 2 in apertures centered on peaks of CO and 
Ha emission, respectively. Unfortunately, they do not 
report the scatter in the £h 2 — Ssfr relation. However, 
we can compare the change of log 10 (tqo — t h 2 ) , a crude 
estimator of the scatter, with changing spatial averaging 
scale using their Fig. 3. This results in a scaling similar 
to the predictions given in our Fig. [4] 

To summarize, we find that for galaxies with Milky- 
Way like ISM conditions most of the scatter in the 
Eh 2 — EgpR relation is a consequence of the time dis- 
creteness of star formation. Systematic variations in the 
SFR tracer conversion factors are only relevant for trac- 
ers with short lifetimes (e.g., Ha) and when observations 
are done on sufficiently small scales (< 400 pc). Varia- 
tions of the H2/CO conversion factor can dominate the 
scatter on kpc scales and above, but are unimportant 
on much smaller scales. Finally, fluctuations in the H2 
abundance play a supporting role, enhancing the scatter 
caused by the discreteness of the star formation process. 

4. DISCUSSION 

How does the star formation model presented in this 
paper relate to alternative interpretations of the scatter 
in the £h 2 — Esfr relation? 

A commonly made suggestion is that the relation 
"breaks down" on small scales ( Momose et al.1 120101 : 



lOnodera et alT20Tbl ISchruba et afl feoiO) because molec- 
ular clouds pass through evolutionary phases in which 
they transform from CO-bright, but star-less, clouds to 
star forming regions with little surrounding molecular 
gas. One of the difficulties with a picture in which 
most of the molecular gas in a galaxy goes through 
a well defined sequence of stages is the following. It 
does not explain why the H2 depletion time of molec- 
ular clouds with embedded young stellar objects is only 
a few 100 Myr (jLada et al.l 12010( 1. while the depletion 
time on galactic s cales is an order of magnitude larger 
(jBigiel et al.l 1201 il l. In fact, this observation, originally 
used as evidence to support long lifetimes of molecular 
clouds (|Zuckerman fc Evanslll974H . implies that the ma- 
jority of the molecular g as in the galaxy has to be in 
a non-star forming state (jElmegreenl l2000t l . possibly in 
either non-star forming clouds that will be dispersed be- 
fore star formation has a chance to begin, or in unbound 
molecular associations. 

The way our model addresses this problem is that the 
observed gas depletion time of a molecular and star form- 
ing region is smaller than the average depletion Td ep by 
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a factor ~ Ai*/AisF (see appendix). Observations of 
individual molecular clouds often derive SFRs based on 
counts of young stellar objects (At* = 1 — 2 Myr) or 
by using tracers of ionizing radiation from massive stars 
(Af» ~ few Myr), e.g., Ha emission or free- free radio 
emission. Consequently, At* < Atgp = 10 Myr and the 
observed gas depletion time in star forming regions is 
shorter than Td cp - This is balanced by the large (for- 
mally infinite) depletion time in molecular regions that 
are not currently star forming. 

It has also been suggested that the scatter in the 
Sjj 2 — Ssfr relation is caused by variations in the ratio 
between molecular gas and dense (> 10 4 cm~ 3 ) molec- 
ular gas (jLada et al.ll2012D . Since our model does not 
explicitly follow dense gas, such variations could be in- 
cluded as an additional stochastic component in the star 
formation model, see Q3A\ However, the tight correlation 
between dense ga s and star formation in many differ- 
ent e nvironments (|Lada et alj|2012t IPapadopoulos etahl 
120H is suggestive of an alternative way. In the context 
of our model the only thing required is to re-interpret 
the word "individual star formation event" as "individual 
dense gas formation event" and to assume that, once gas 
becomes very dense, star formation is inevitable and will 
proceed on the ~ free-fall time of the respective dense gas 
clump. This modification of our model does not specify 
the physical mechanism for the sudden increase in gas 
density, b ut several plausible options exist, e.g., cloud 
collisions (jTasker fc Tanll2009ft . 

This re-interpretation of our star formation model ac- 
counts, by construction, for the observed tight correla- 
tion between dense gas and star formation rate. Fur- 
thermore, the ratio between dense (i.e., star forming) 
molecular gas and all molecular gas in a given region 
will vary depending on how many "dense gas formation 
events" have occurred in the region. In this model both 
the scatter in the £h 2 — ^SFR relation and the scatter 
in the mass ratio between dense and not-so-dense gas on 
small scales are a consequence of the discrete formation 
of dense gas clumps out of molecular gas. 

Further potential contributors to the scatter include 
the incomplete sampling of the IMF and the drifting of 
stars out of their parent molecular clouds. These scat- 
ter sources are unlikely to be relevant given the rela- 
tively large spatial scales (> 100 pc) and SFR surface 
densities (> 10~ 3 M yr _1 kpc -2 ) in our study (see 
lOnodera et alll2010h . 

We conclude that in the context of our star formation 
model there is no "break down" of the Eh 2 — ^sfr scaling 
relation. Instead, the proper interpretation is that the 
discrete and stochastic nature of star formation becomes 
evident as observations probe smaller and smaller scales. 

5. SUMMARY AND CONCLUSIONS 

In this paper we studied the slope and the scatter of 
the Sh 2 — ^sfr relation using cosmological galaxy forma- 
tion simulations coupled with models for star formation, 
H2 chemistry, and CO emission. Our focus is especially 
on the role of the CO/H2 conversion factor. We found 
that X-factor variations with surface density can result 
in significant biases of the measured slope. In particular, 
at high spatial resolution (few 100 pc or better) and suf- 
ficiently low surface densities (Sh 2 < 20M Q pc -2 ) the 
slope inferred from CO observations is shallower than 



the actual slope. In contrast, the inferred slope becomes 
steeper than the true slope if galaxies with high H2 sur- 
face densities (above 100 Mq pc -2 ) are included in the 
sample, providing an possible explanation for the slightly 
super-linear slope of the S h q — £sfr relation se en at high 
gas surface densities (e.g.. iGenzel et al1l2010f ). Yet, we 
also showed that measurements at > 500 pc resolution 
over a surface density range often studied in samples of 
nearby galaxies, 10M Q pc~ 2 < Er 2 < 100Af©pc -2 , are 
essentially unbiased. 

Variations in the X-factor contribute to the scatter in 
the Eh 2 — Ssfr relation (of the order of ~ 0.1 — 0.2 dex), 
dominating over many other scatter sources when the 
spatial resolution of the survey is ~ kpc or larger. This 
even holds if there are no significant spatial variations 
of the dust-to-gas ratio or the interstellar radiation field. 
Such variations are expected in a heterogeneous sample 
of galaxies, leading to additional galaxy-to-galaxy scat- 
ter with an amount that depen ds on the properties o f the 
particular galaxy sample (e.g.. iSchruba et al]|201lD . On 
sub- kpc scales, however, spatial variations in the CO/H2 
conversion factor contribute little to the overall scatter 
(assuming a fixed dust-to-gas ratio). On such scales 
much of the scatter is a consequence of the fact that 
the measured, time-averaged SFRs differ from the SFRs 
that are expected based on the present amount of H 2 . 

We demonstrated that the scatter in the £h 2 — Ssfr 
relation (on scales of ~ 100 pc and larger) is primarily 
a consequence of the discreteness of the star formation 
process. The luminosity evolution of SFR tracers can 
become relevant for tracer with short lifetimes, e.g., Ha, 
on small scales (less then a few 100 pc). For FUV-based 
SFRs, however, the scatter does not differ significantly 
from the scatter based on a hypothetical tracer with a 
plain 20 Myr lifetime. The differences in the timescales 
between SFR tracers (at least a few Myr, up to 100 Myr) 
and that of H2 masses (essentially instantaneous mea- 
surements) does lead to some scatter, but it is typically 
dominated by scatter that results from X-factor varia- 
tions (on super-kpc scales) and by scatter due to the 
discreteness of star formation (on sub-kpc scales). 

The predictions made in the paper suggest a number 
of observational tests that could be used to constrain the 
presented numerical models. Fortunately, most of the 
more obvious tests, e.g., looking for a change in slope of 
the Eh 2 — Ssfr relation at very low and very high surface 
densities (Fig 03, Fig [3]), should be feasible with future 
ALMA observation. A clear test of the predictions of the 
Poisson star formation model should be possible with a 
systematic, observational study of how the scatter in the 
Eh 2 — Ssfr relation scales with spatial scale and how it 
depends on the lifetimes of star formation tracers. Such 
a study would allow to constrain the discreteness of star 
formation and, more generally, would be a crucial guide 
for the development of the theoretical underpinnings of 
star formation in a galactic context. 
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APPENDIX 

In this appendix we give a more formal definition of the Poisson star formation model introduced in £ !2.2l and describe 
its implementation in the ART code. The model assumes that the number of individual star formation events N At in 
a given time interval At is a Poisson random variable with a mean and variance A = ■ Here, Aigp is the average 
time interval between two star formation events. The continuum, non-stochastic limit corresponds to AtsF 0. 

The ensemble average SFR density in a given region is proportional to the H2 density, i.e., 

</>♦>(*) = —, W 

7"dcp 

where the gas depletion time Td ep is assumed to be a constant. The ensemble average SFR in a given fixed volume Z| F 
is simply S = (p*)^g F . We denote the time average of S over some interval At as S A t. 
The actual SFR, Sa*, that occurred in the volume Z| F during the time interval At is a random variable and defined 

as 

S At = N At ^S Au (2) 

which ensures that the ensemble average of S A t equals S A t- This equations shows that the stellar mass formed during 
an individual star formation event is simply given by S A tAtsF, i-e., it scales linearly with the average SFR. 

The time interval At may be some small, arbitrarily chosen time step. We are typically interested in measuring the 
mean and variance of the SFR (or of its logarithm) over some physical time interval At* 3> At, e.g., over the lifetimes 
of star formation tracers such as Ha or FUV emission. The actual SFR over At* is given by 

Sau = -^r -Sam = -^p- N M,i S A t,i, (3) 

* i * i 

with indices i running over the At* /At time intervals of length At. 

It is worthwhile to have a closer look at ((3]). First, the expression on the right hand side implies that Sau does not 
depend on the time step At provided At is sufficiently short compared with the typical time over which S fluctuates. 
Moments of 5a*, and log 10 5a*, depend on the timescales AtsF and At* only via the ratio At*/Atgp. The equation 
also allows us to estimate the H 2 depletion time pu 2 ^sf/Sau that an observer would measure in a star forming region 
(TVa*. = N A t.i > !)■ If we assume that the H2 density in the region remains constant during At*, then the 
depletion time that the observer would infer is Tdcp(At t / Atg-p)/N Atr , which is smaller than Td ep if At* < AtgF- 

Equation ([3]) further shows that Sau depends in a non-linear way on both the time evolution of the H2 mass in a 
given region and the number of individual star formation events. Hence, we expect that the scatter in the Sr 2 — ^SFR 
relation that arises from the discreteness of star formation depends in a non-linear way on both the time variations 
in the H 2 mass and on Poisson shot noise in the number of individual star formation events. To be quantitative let 
us approximate the scatter in the £h 2 — Esfr relation with ai og s At > the standard deviation of the logarithm of 
Sau > 0- Let us further assume that Sam fluctuates over time At* as a log-normal random variable, i.e., Sam 00 e<yX 1 
where A is a Gaussian random variable with mean and variance 1. The choice At* = 20 Myr, AtsF = 10 Myr and 
a = 1.4 results in cri ogl0 s At , ~ 0.60. With the same set of parameters but AtsF — > we find cr logio g Atm ~ 0.07. If we 
assume a — 0, i.e., 6>am constant, we obtain ci ogio s& t , ~ 0.24. Hence, the variance of log 10 Sau is n °t simply the 
sum of the variances caused by Poisson noise and H 2 fluctuations, respectively, but is determined to a large extent by 
their covariance. Using this simple test setup we also find that the scaling of <xiog 10 SAt, with At*/AtsF depends on 
the value of a, although, for sufficiently large At*/AtsF, the scatter decreases with increasing At*/AtsF- 

Our simulations adopt the Poisson star formation model in the following way. Every At = 10 5 yr the code computes 
the ensemble average SFR S in each resolution element (1$f ~ 100 pc, see Table [T]) based on the present H 2 mass. 
Then, a random realization of N At is drawn from a Poisson distribution with the mean At/AtsF and, if N At > 0, the 
code creates a stellar particle of mass N At <SAtsF ■ 
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